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ABSTRACT 

We present detailed clustering measurements for a flux limited sample of 14, 000 
quasars extracted from the 2dF QSO Redshift Survey (2QZ) in the redshift range 
0.8 < z < 2.1. After splitting the sample into three redshift bins and each of them into 
six luminosity intervals, we estimate the quasar projected auto and cross-correlation 
functions at a given redshift for separations 3 < r//i^^Mpc < 20. Fitting the data with 
a biased CDM model and using a frequentist analysis (the F-test), we find that models 
with luminosity dependent clustering are statistically favoured at the 95 per cent 
confidence level for z > 1.3. On the other hand, a number of tests based on information 
theory and Bayesian statistics show only marginal evidence for luminosity dependent 
clustering. Anyway, the quality of the data is not good enough to accurately quantify 
how quasar biasing depends on luminosity. We critically discuss the limitations of our 
dataset and show that a much larger sample is needed to rule out current models 
for luminosity segregation. Studying the evolution of the clustering amplitude with 
redshift, we detect an increase of the quasar correlation length with lookback time at 
the 99.3 per cent confidence level. Adopting the concordance cosmological model, we 
discuss the evolution of quasar biasing with cosmic epoch and show that quasars are 
typically hosted by dark matter haloes with mass ^ 10^'^AfQ. 

Key words: galaxies: active - galaxies: clustering - quasars: general - cosmology: 
theory - large-scale structure - cosmology: observations 



1 INTRODUCTION 

It is widely believed that quasars are powered by accretion 
onto supermassive black holes. However, a detailed under- 
standing of the physical processes leading to quasar activity 
(and their connection with galaxy formation) is still lacking. 

Simple semi-analytic models associate quasars with 
galaxy major mergers and assume a tight relation between 
their instantaneous luminosity and the mass of the central 
black-hole, Mbh (Kauffmann & Haehnelt 2000; Wyithe & 
Loeb 2003; Volonteri et al. 2003). The fraction of gas ac- 
creted onto the black hole during each merger is chosen to 
match the observed relation between the velocity dispersion 
of the bulge and Mbh (Ferrarese & Merritt 2000) . This ends 
up producing a correlation between the quasar luminosity 
and the mass of the host dark-matter halo. Since the clus- 
tering properties of dark-matter halos strongly depend on 
their mass, the quasar clustering amplitude is thus expected 
to sensibly depend on luminosity. 

Recent numerical simulations of galaxy mergers includ- 
ing black-hole accretion and feedback have cast some doubts 
on this picture (Springel et al. 2005). These numerical ex- 



The Scottish Universities Physics Alliance 



periments suggest that a given black-hole produces quasar 
activity with a wide range of luminosities (Hopkins et al. 
2005). During its active phase, the black- hole is most likely 
observed as a relatively low-luminosity quasar with a small 
Eddington ratio. For a short period of time, however, its 
emission reaches its peak value (close to the Eddington lumi- 
nosity) which is indeed proportional to the mass of the pow- 
ering black-hole. Based on these models, Lidz et al. (2006) 
conclude that quasar clustering should depend only weakly 
on luminosity. 

From the observational point of view, only recently 
quasar samples have grown big enough (in terms of number 
of objects) to attempt the study of the clustering amplitude 
as a function of luminosity. By analyzing the galaxy-quasar 
cross-correlation at 1.8 2: < 3.5 Adelberger & Steidel 
(2005) found no evidence for luminosity-dependent cluster- 
ing. They used 79 quasars spanning 4.4 orders of magnitude 
in absolute luminosity which have been divided into 2 lumi- 
nosity bins. Larger samples are obviously required to confirm 
this result. 

Croom et al. (2002, 2005) studied the redshift-space 
clustering amplitude of 2dF quasars as a function of their 
apparent magnitude. Even though these authors initially 
found weak evidence for brighter QSOs being more strongly 
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clustered, their most recent analysis shows no indication of 
luminosity-dependent clustering. These two studies, how- 
ever, do not address the issue of luminosity dependent clus- 
tering. In fact, they consider magnitude- limited samples 
within a broad redshift range (0.3 < 2 < 2.9 corresponding 
to ~ 8 Gyr in the currently favoured cosmological model) 
and totally ignore any changes in the clustering signal as 
a function of cosmic epoch. Moreover, their redshift-space 
analysis complicates the interpretation of the clustering am- 
plitudes, as the effect of non-linear peculiar velocities could 
also depend on luminosity and redshift. For these reasons, 
their null result does not necessarily imply that quasar clus- 
tering is independent of luminosity. 

In this paper, which is a follow up analysis of Porciani, 
Magiiocchetti & Norberg (2004, PMN04 hereafter), we study 
the clustering properties of ~ 14, 000 quasars extracted from 
the complete 2dF QSO Redshift Survey (Groom et al. 2004). 
Our goal is to accurately measure the real-space clustering 
amplitude of 2dF quasars as a function of redshift and abso- 
lute luminosity. To do this, we first split our quasar sample 
into three redshift bins and each subsample into six luminos- 
ity intervals. We then compute, in each redshift range, the 
associated projected auto- and cross-correlation functions, 
which, by construction, are not affected by redshift space 
distortions. 

Using the largest quasar sample presently available, our 
results reveal a statistically significant evolution of the clus- 
tering length with redshift (as found already in PMN04) but 
only weak evidence for luminosity dependent quasar cluster- 
ing. 

The layout of the paper is as follows. In Section |21 we 
describe our quasar sample and how we split it both in red- 
shift and in luminosity. In Section |21 we measure the pro- 
jected auto- and cross-correlation functions of the different 
subsamples. The issue of luminosity dependent clustering is 
addressed in Section |1| Using a Monte Carlo Markov Chain, 
we estimate the quasar clustering amplitude as a function of 
redshift and luminosity. A number of robust tests are then 
used to evaluate the statistical significance of the measured 
luminosity dependence in a given redshift range. A critical 
discussion of the limitations of our data and of possible fu- 
ture strategies is also presented here. In Section |3| we focus 
on the pure redshift dependence of the quasar clustering am- 
plitude and we provide fitting functions for the evolution of 
the quasar correlation length and bias with cosmic epoch. 
All our results are summarized in Sectional 

Throughout this paper we assume a "concordance" cos- 
mological model with mass density parameter fio — 0.3 
(with a baryonic contribution Qb = 0.049), vacuum en- 
ergy density parameter flA = 0.7 and present-day value of 
the Hubble constant Hg = 100 /ikms"'^ Mpc with h = 0.7. 
We also adopt a cold dark matter (GDM) power spectrum 
with primordial spectral index n = 1 and with normaliza- 
tion fixed by as, the rms linear density fiuctuation within a 
sphere with radius of 8 Mpc. 



2 QUASAR SAMPLE DEFINITION 

The 2dF QSO Redshift Survey (2QZ) includes 23,338 
quasars which span a wide redshift range (0.3 ^ z ^ 2.9) 
and are spread over 721.6 deg^ on the sky (see Groom et al. 



Table 1. Main properties of our data sets. The subscripts min, 
max and med respectively denote the minimum, maximum and 
median value of a variable. 
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2004). In order to minimise systematic effects, we restrict our 
analysis to regions with spectroscopic (photometric) com- 
pleteness larger than 70 (90) per cent, which limits the red- 
shift range to 0.5 < z < 2.1. Only quasars brighter than 
M(,, — Slogj^Q h — —21.7 are considered, which ensures the 
exclusion of quasars for which the contribution from the host 
galaxy may have led to a mis-identification of the source. 

In order to make a physically motivated analysis and 
isolate evolutionary effects, we subdivide our sample into 
three redshift bins. As in PMN04, we require a similar num- 
ber of quasars to lie in each redshift bin and that each sub- 
sample covers a not too different interval of cosmic time. 
To better satisfy these criteria (see PMN04 fur further de- 
tails), we impose an additional redshift cut so as to keep 
only quasars within 0.8 < z < 2.1. With this selection, we 
end up with nearly 14,000 2QZ quasars that we split into 
the redshift intervals: 0.8 < z < 1.3, 1.3 < z < 1.7 and 
1.7 < z < 2.1 (containing each between ~ 4300 and ~ 4900 
quasars) . 

To study the dependence of the quasar clustering ampli- 
tude on luminosity, we further divide each sub-sample into 
six complementary sets based on quasar absolute luminos- 
ity: (F25, B75), (F50, B50) and (F75, B25), where Fx and 
By correspond to the x per cent faintest quasars and the y 
per cent brightest quasars (in terms of their absolute mag- 
nitude). Table lists the main properties of each sample, 
including the number of quasars (col. 4), the median red- 
shift (col. 5), the absolute magnitude range (cols. 6 & 8) 
and the sample median absolute magnitude (col. 7). 
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Figure 1. Top-left panel: The projected auto-corrclation function (top two rows) and cross-correlation function (bottom two rows) for 
the low-redshift quasar samples (0.8 < z < 1.3). Only measurements obtained with more than 40 quasar pairs are plotted (filled squares 
with errorbars). Bins with less pairs are omitted and not considered for model fitting. For each luminosity sub-sample, the shaded band 
show the central 68 per cent range of best-fitting correlation functions (assuming that each luminosity interval is associated with a 
difli'erent bias parameter). The mean value of the marginalized bias distribution is represented with a solid line. For reference, we also 
plot the best- fitting correlation function for the entire low-redshift sample derived in PNM04 (dotted line). Top-right and bottom-left 
panels: Same as the top-left panel, but for the medium- (1.3 < z < 1.7) and high- (1.7 < z < 2.1) redshift quasar samples respectively. 
Bottom-right panel: As in the top-right panel but using a three-parameter model where only the B25 and F25 sub-samples have different 
bias parameters from all the others. 



3 QUASARS CLUSTERING ANALYSIS 

3.1 Estimating the correlation functions 

The simplest statistic which can be used to quantify clus- 
tering in tiie observed quasar distribution is the 2-point cor- 
relation function in redshift space, ^''(rx,7r). To measure 
this quantity we first build a catalogue of unclustered points 
which has the same angular and radial selection function as 
the data. The radial selection function is obtained by heavily 



smoothing the observed quasar comoving distance distribu- 
tion, M{r). The quasar auto-correlation function is then es- 
timated by comparing the probability distribution of quasar 
and random pairs on a two-dimensional grid of separations 
(rx, 71"). We use both the Landy-Szalay estimator (Landy & 
Szalay 1993) and the Hamilton estimator (Hamilton 1993): 

_DD-2DR + RR _ DD ■ RR 

- RR ' " {DRY ~ ^ ' 
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where DD, DR and RR are the suitably normahsed numbers 
of weighted data-data, data-random and random-random 
pairs in each bin/ As expected, the two estimators give 
comparable answers within the errors. In what follows we 
only present results obtained with the Landy-Szalay estima- 
tor. Similarly, we use the appropriately symmetrized Landy- 
Szalay estimator to compute the quasar cross-correlation 
function: 

q _ DjDj — {DjRj + DjRj) + RjRj , , 

tXittj 

where DiDj, DiRj, DjRi and RiRj are the suitably nor- 
malised numbers of weighted data and random pairs from 
samples i and j. 

To avoid redshift space distortions, and determine the 
quasar clustering amplitude in real space, one can then use 
the 'projected correlation function' which is obtained by in- 
tegrating ^''('■-Li^) in the tt direction: 

^ = ^ re''(rx,^)d.. (3) 

With the current quasar sample, we find that a reliable 
measure of ^'^{r±,n) is only achievable on scales vr < 
50 Mpc. In order to avoid the measured signal to be 
dominated by noise, we limit the integration to an upper 
limit, TTmax. As explained in PMN04, TTmax = 45 h'^ Mpc 
is an adequate choice, and we have checked that our results 
are not sensitive to the precise value adopted. 

As in PMN04, we use a blockwise bootstrap resam- 
pling technique to estimate the uncertainties of our clus- 
tering measurements. For each redshift interval, we divide 
each of the two 2QZ areas into 28 equal-volume regions, con- 
taining, within a factor of two, the same number of quasars 
(the minimum size of a subvolume is nearly lSOh~^ Mpc 
for the low-redshift sample, 250^"^ Mpc for the medium- 
redshift sample and 220h~^ Mpc for the high-redshift; since 
we measure correlations on scales smaller than 20h~^ Mpc 
our error estimates should hardly be affected by edge-effects 
and/or missing large-scale structure). ^ We build hundreds 
of bootstrap-samples (see tl3.31 , each of them composed by 
56 sub-samples (28 for each 2QZ region) randomly drawn 
(allowing repetitions) from the set described above. We mea- 
sure the projected auto- and cross-correlation functions for 
all the bootstrap samples. For each r±, we identify the rms 
variation of H'' over the bootstrap-samples with the la error 
for the projected auto- and cross-correlation function. Note 
the our bootstrap uncertainties account for both sampling 
and estimation errors. 

Our method for estimating errors relies on the fact that 
our dataset is statistically representative of the quasar dis- 
tribution in the Universe. ^ However, this cannot be true for 

^ Note that, in this case, there is no need to use the standard J3 
(minimum variance) weighting scheme since the mean density of 
quasars, Uq, is so low that 1 -|- 4 tt J3 ?iq ~ 1 for any reasonable 
quasar clustering amplitude. 

^ We have also checked that doubling the size of the subvolumes 
does not significantly alter the diagonal errors. 

Our blockwise bootstrap errors are in good agreement with 
those obtained using mock galaxy catalogues extracted from high- 
resolution numerical simulations. Discrepancies between the two 
types of errors are of order 20 per cent or less on the scales used 
in this paper. 



bins of spatial separations which contain just a few quasar 
pairs. Therefore, in what follows, we ignore clustering results 
obtained with less than 40 quasar pairs. Depending on the 
sample, this corresponds to < 2 — 6 h^^ Mpc. 

Our results are presented in Fig. (filled squares with 
errorbars). All the correlations in a given redshift bin have 
similar amplitudes and are rapidly decreasing with r±. As 
expected, errorbars are smaller for the largest quasar sam- 
ples. 

3.2 Modelling the data 

Given the relatively large errorbars, it is hard to spot any 
luminosity dependence of quasar clustering by simply look- 
ing at the correlations in Fig. For this reason, we use a 
reference model to fit the data. This is obtained by assuming 
that each sample is characterized by a linear bias parame- 
ter which does not depend on spatial separation within the 
range of scales we analyse. In other words, we assume that 
the quasar 2-point auto- and cross-correlations scale as 

Q{r) = b,b,ar), (4) 

where i and j axe the labels of 2 quasar sub-samples within 
the same redshift bin, r denotes the comoving separation be- 
tween quasar pairs, and bi is the bias parameter of the i-th 
sub-sample with respect to the mass autocorrelation func- 
tion, 5(r), computed as in Peacock & Dodds (1996, nearly 
undistinguishable results are found using the method by 
Smith et al. 2003) assuming a value for erg. * Within the 
framework of halo models (see PMN04 for a direct applica- 
tion to quasars), the assumption that the amplitude of the 
cross-correlation function between haloes of different masses 
scales as the geometric mean of the correlations of the in- 
dividual haloes is very natural in the two-halo regime. We 
have checked against numerical simulations that, for dark 
matter haloes, the assumption holds over the redshift range 
and scales considered here, i.e. 3 < r < 20 Mpc. The 
corresponding projected correlation functions are obtained 
through this simple integral relation: 

= 2 dr. (5) 

3.3 Fitting the data with correlated errors 

We use a minimum least-squares method (corresponding to 
a maximum likelihood method in the case of Gaussian er- 
rors) to determine the bias parameters that best describe the 
clustering data. In each redshift range, we sample the six- 
dimensional parameter space (one bias parameter per lumi- 
nosity range) using a Markov Chain Monte Carlo method. A 
principal component analysis (see e.g. Porciani & Giavalisco 
2002) is used here to deal with correlated errorbars. The 

* Note that samples with different luminosities have slightly dif- 
ferent median redshifts. However the variation of the mass corre- 
lation function among them is of the order of a few per cent, much 
smaller than the uncertainty in the quasar clustering amplitude. 
This implies that we can safely use the same ^(r) (evaluated at 
the median redshift of the redshift bin) for all the luminosity sub- 
samples. We have tested that this does not influence our results. 
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principal components of the errors are computed by diago- 
nalizing the covariance matrix obtained by resamphng the 
data with the bootstrap method. The objective function (the 
usual statistic) is then obtained by only considering the 
most significant principal components (i.e. those contribut- 
ing the largest fraction of the total variance). ^ Unfortu- 
nately there is no unique objective way of deciding how 
many principal components should be considered for a given 
dataset. Considering too few components (and thus discard- 
ing information contained in the data), one obtains too good 
values for the objective function (i.e. Xmin/dof <^ 1, where 
dof indicates the number of degrees of freedom) that cor- 
respond to unrealistically large uncertainties for the fitted 
parameters. On the other hand, considering too many com- 
ponents (and thus, most likely, introducing noise), one of- 
ten obtains bad fits (i.e. Xmin/dof ^ 1) that correspond 
to unrealistically small errors for the bias parameters. In 
factor analysis, a number of empirical methods (e.g. Kaiser 
criterion, scree test) are often employed to select the num- 
ber of significant components. The robustness of these tech- 
niques for model fitting is, however, rather weak. We thus 
decided to select the number of components to account for 
a fixed fraction (95 per cent) of the bootstrap variance. We 
have checked that this compression guarantees a detailed 
reconstruction of the model-data residuals and simultane- 
ously avoid the to be dominated by deviations along 
the components corresponding to tiny eigenvalues. More- 
over, the reduced of the best-fitting models obtained this 
way is of order unity in all cases, as expected theoretically. 
In what follows, we thus use the symbol X95 to denote the 
X^ function computed by considering the first principal com- 
ponents that, in total, account for 95 per cent of the vari- 
ance. This corresponds to 19 principal components (out of 
42 datapoints) for the high-redshift sample, 21 (out of 45) 
for the median-redshift sample and 22 (out of 47) for the 
low-redshift sample. 

Since the bootstrap technique is very time consuming, 
a careful choice of the number of resamplings is required. 
For this reason, we performed a number of Monte Carlo 
simulations checking for the stability of eigenvalues, eigen- 
vectors and X95 estimates. In practice, we first bootstrapped 
our data and computed the corresponding covariance matrix 
Cij. We then built a large number of realizations of a Gaus- 
sian process with covariance dj. Finally, we estimated the 
covariance matrix of the Gaussian process from a finite num- 
ber of realizations, A^, and studied its dependence on A'^. We 
found that a few hundred bootstrap resamphngs are needed 
for a robust estimation of the X95 function. To overcome the 
Gaussian hypothesis, we also studied the convergence of the 
covariance matrix obtained by directly bootstrapping our 
data. We found that estimates of X95 converge (i.e. present 
neghgible scatter) when nearly 250 bootstrap resamplings 
are used. This has been obtained by using 500 resamplings 
of the median redshift sample. To be on the safe side, we 

^ The covariance matrix obtained from bootstrap resampling is 
only an estimate of the true one and contains an intrinsic uncer- 
tainty. The errors in its component propagate in the calculation 
of its eigenvalues and eigenvectors. It is therefore recommended 
to consider only the principal components corresponding to the 
largest eigenvectors in the fitting procedure (see §4.2 in Porciani 
& Giavalisco (2002) for further details). 



thus decided to use at least 350 bootstrap resamplings for 
each redshift bin. We note that this implies over 12,000 cor- 
relation function estimates, with each containing at least 
50,000 random points. 

As an additional test of the robustness of our results, 
we checked, a posteriori, how much our best-fitting models 
depend on the number of adopted components, A^pca. This is 
discussed in detail in Section^Jwhere we present our results. 



4 LUMINOSITY DEPENDENT CLUSTERING 

In this section, we present the results obtained by fitting 
the auto- and cross-correlation functions presented in Fig. Q 
with the model given in equation Our aim is to quantify 
how the bias parameter of optically bright quasars depends 
on redshift and luminosity. 

The best-fitting functions obtained with the Markov 
Chain Monte Carlo method are overplotted to the data in 
Fig.0(shaded regions). Note that, for all redshift bins, they 
accurately describe the quasar auto- and cross-correlations 
in the whole range of separations under analysis. In the left 
panel of Fig. |21 we plot cTg (i.e. bag) as a function of ab- 
solute magnitude. As expected, the bias parameter of all 
quasars lying in a given redshift bin (measured by PMN04 
and represented with an open symbol) lies in between the 
values found here. Taken at face value, the data show some 
evidence for luminosity segregation; however, samples at dif- 
ferent redshifts show different trends. The brightest quasars 
in the high-redshift bin seem to be more strongly clustered 
than the others. In the medium-redshift bin, the bias seems 
to follow a [/-shape: the faintest and the brightest quasars 
in the set are typically more strongly clustered than MJ^ 
quasars. On the other hand, the low-redshift bin does not 
show any particular trend: the bias keeps constant with lu- 
minosity. 

Note, however, that correlations between errorbars in 
Fig. m might create spurious trends with luminosity. More- 
over, all the different sub-samples in a given redshift bin have 
consistent bias parameters at the 2a level. Given these con- 
cerns, we want to use a robust statistical test to investigate 
whether our data show any sign of luminosity dependent 
clustering. 

4.1 Does quasar clustering depend on luminosity? 

If there is no luminosity segregation, all our data should be 
described by one single bias parameter (i.e. all the bi should 
assume the same identical value) . If, instead, some luminos- 
ity bins show a statistically significant deviation from the 
overall clustering amplitude, then a description with a imm- 
ber of different bias parameters should be preferred. There- 
fore, in this section we address the issue of luminosity de- 
pendent clustering by answering the following question: how 
many bias parameters are required to adequately describe 
the quasar auto- and cross-correlations in a given redshift 
bin? In other words, we want to understand how many dif- 
ferent bias parameters can be reliably measured in a statis- 
tically significant way. This is a classical problem of model 
selection where we want to find the proper tradeoff between 
goodness of fit (in the x^ sense) and complexity (in terms of 
number of free parameters). We use five different methods: 
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Figure 2. Left panel: cr^, i.e. bag, as function of the median absolute magnitude for the low (left), medium (centre) and high (right) 
redshift bins. For each sample, the errorbars indicate the central 68.3 (solid) and 95.4 (dotted) credibility intervals of the marginalized 
bias probability distribution. On the upper axis, we indicate the magnitude difference with respect to MJ^(z), where the later is given 
by Eq. 13 of Croom et al. (2004). The open symbol show the result obtained by fitting the correlation function of all the quasars lying 
in this redshift bin (see PMN04). Right panel: Same as the left panel but for a three-parameter fit (see text). 



the F-test, the Akaike information criterion (AIC), the Bayes 
factors (BF), the Bayesian information criterion (BIC) and 
the deviance information criterion (DIG, see the Appendix 
for a brief review). We consider three different models. In the 
first one, there is no clustering segregation with luminosity 
and all the luminosity bins at a given redshift are associated 
with the same bias parameter. In the second model, we use 
three bias parameters (one each for the brightest, B25, and 
faintest samples, F25, and one for the remaining quasars: 
hereafter those samples are referred to as B, M and F, for 
bright, medium and faint). In the third model, each luminos- 
ity bin has its own bias parameter, as already described in 
the previous sections. In Fig. |21 we show that the results of 
the three- and six-parameters models are qualitatively and 
quantitatively similar. 

Before comparing the different models, it is worth notic- 
ing that, based on the statistic, a single parameter fit 
gives an acceptable description of the data for all redshift 
bins. Of course, using additional parameters reduces the 
minimum of the best-fitting models. We want now to 
understand whether this Xmin reduction is statistically sig- 
nificant. 

Performing the _F-test and using the 95 per cent con- 
fidence level as a threshold to prefer a model with respect 
to another, we find that the low-redshift sample is best- 
fitted by a one-parameter model, while the medium- and 
high-redshift samples are best described by a six-parameter 
model. Similarly, all the Bayesian and information-theory 
based tests mentioned above clearly indicate that the low- 
redshift sample is best described by a one-parameter model. 
On the other hand, for the medium- and high-redshift sam- 
ples, the situation is more confused. Because of the different 



Table 2 . Number of parameters of the preferred models according 
to various model-selection criteria. When the evidence for the best 
model is not strong, we list all the acceptable models (in order 
of decreasing evidence). We use a 95 per cent confidence level for 
the F-test, AIC and DIG while we discard models with substantial 
evidence according to the Kass-Raftery criterion (see Appendix) 
for BF and BIC. 



Criterion O.S < z < 1.3 1.3 < z < 1.7 1.7 < z < 2.1 



F-test 16 6 

AIC 1-3 1-3 1-3 

BF 1-3-6 6-3-1 6-3-1 

BIC 1 1-3-6 1-3-6 

Die 1-3-6 6-3-1 6-3-1 



penalties for model-complexity, depending on the adopted 
test, either the six- or the one-parameter model is the pre- 
ferred one. However, with the exception of the six-parameter 
fit with the Akaike criterion, no model can be rejected with 
high- confidence (see Table |5| for a summary of the results). 
Therefore, we conclude that, while there is no evidence of 
clustering segregation with luminosity in the low-redshift 
sample, we find marginal evidence for it in the medium- and 
high-redshift samples. 

Following a suggestion of the referee, we have also tried 
to combine data from the different redshift bins by us- 
ing a simple parameterization of the luminosity and red- 
shift dependent bias. For simplicity, we adopted a lin- 
ear relation with quasar luminosity b{Mi,j,z) = bo{z) + 
biW ''■'"^^^''.j jjj ^jjjg case, the 68.3 per cent con- 

fidence interval for the luminosity-dependent parameter is 
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bi — O-Slp y. This confirms tiiat the evidence for cluster- 
ing segregation with luminosity is marginal in our sample. 
Similar conclusions are drawn adopting quadratic or cubic 
relations for the luminosity dependence of the bias parame- 
ter. In all cases, the best-fitting model favours a larger bias 
parameter for the brightest quasars in our sample. However, 
the statistical significance of the result is low. 



4.2 Limitations of the PCA method 

A number of systematic efi'ects can alter the model-fitting 
procedure (and thus the statistical significance of the re- 
sults). In particular, for strongly correlated data, the choice 
of the number of principal components used in the fit plays 
a delicate role. We fixed A'pca assuming that our bootstrap 
errors are accurate estimates of the real uncertainties and 
using the test. Had we decided to include a few more com- 
ponents to account for 99 per cent of the variance (which is 
equivalent to consider 7 additional components for each red- 
shift sample and still gives acceptable values of the reduced 
X^), the evidence for clustering segregation would have been 
stronger. In Fig. |31 we compare the solutions for the three- 
parameter models obtained by minimizing xis (left) and 
X99 (right). What is shown here are the (marginal) joint 
probability distributions of pairs of bias parameters (for the 
faint, medium and bright samples). Clustering segregation 
with luminosity is suggested by the data whenever the con- 
tours of the joint distributions lie away from the diagonal 
line. For both xis and xlgi this never happens for the low- 
redshift sample (the bottom-left set of curves) and the one- 
parameter solution has to be preferred. At best, clustering 
segregation with luminosity is detected at the ~ 2a level (for 
X95) and at ~ 3a level (for X99) in the high- and medium- 
redshift samples. Note that, even though a 3a detection is 
still consistent with pure statistical fiuctuations, when X99 
is used, all the tests for model selection strongly prefer ei- 
ther the three- or the six-parameter solution with respect to 
the one-parameter fit. In other words, the results presented 
in the previous section have to be considered conservative. 
If one decides to use the full covariance matrix, one will 
infer the presence of a statistically significant clustering seg- 
regation with luminosity in the medium- and high-redshift 
samples. Hence, the uncertainty in the number of physical 
components is the main limitation of the principal compo- 
nent technique that we used to account for correlated er- 
rorbars. Bayesian techniques (e.g. Minka 2000) suggest that 
the components which contribute the last few per cent of 
the bootstrap variance are most likely dominated by noise. 
Therefore, we are confident that the results we presented in 
section 4.1 are optimal and realistic. Larger datasets with 
smaller intrinsic uncertainties are thus needed to improve 
the significance of our results. 

In order to cross-check our results with a method 
affected by different systematics, we have computed the 
marked correlation function of the quasars (e.g. Sheth, Con- 
nolly & Skibba 2006 and references therein) using their blue 
luminosity as the mark. Results for the medium redshift 
sample are shown in Fig. 0| Here we define the marked cor- 
relation function (filled squares) as the ratio of the projected 
quasar correlation function weighted by the quasar luminos- 
ity (numerator) and the usual projected quasar correlation 
function (denominator). By construction, this ratio is not 




Figure 4. Marked correlation function of the medium redshift 
sample obtained using the blue luminosity of the quasars as the 
mark. The empty squares (that have been slightly displaced to 
the right) show the mean obtained by averaging over one hun- 
dred realizations of the same measurement after randomizing the 
marks. Error bars show the standard deviation around this mean. 
If one includes bootstrap errors into the total error budget (not 
shown here), the errors on the marked correlation function typi- 
cally increase by 20-30 per cent. 



affected by redshift space distortions. The open symbols cor- 
respond to the mean marked correlation function obtained 
by averaging over one hundred realisations with randomized 
marks. To assess the significance of the marked statistics, 
one usually compares the signal with the standard devia- 
tion around the mean of the randomized realisations. From 
Fig. 4, we see that, on average, the closest quasar pairs have 
higher luminosities, thus suggesting the presence, on those 
scales, of some clustering segregation. However, given the 
size of the errorbars, the statistical significance of this trend 
is rather small. We note that the errors used here only reflect 
the uncertainty arising from the distribution of the marks 
and do not include any uncertainty in the actual estimate 
of the correlation function. Considering the corresponding 
bootstrap errors (which, for instance, allow for sample vari- 
ance), the uncertainties on the marked correlation function 
typically increase by 20-30 per cent. This makes the detec- 
tion of clustering segregation even more uncertain. In sum- 
mary, the analysis of the marked statistics reaches exactly 
the same conclusions as our main work: even though there 
is some evidence for luminosity-dependent clustering in our 
high-redshift data, our sample is too small to provide a ro- 
bust detection. Similar results are obtained for the other 
redshift bins. 



4.3 Discussion 

When interpreting our results, we have to consider the lim- 
itations of both the dataset and the method used to extract 
the bias parameters. First, the depth of the quasar sample 
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Figure 3. Left panel: Results obtained by fitting the data with a three-parameter model (see text) where 63 1 a^nd bp are the bias 
parameters for the bright, medium and faint sub-samples (at fixed redshift), respectively. The contours indicate the 68.3 and 95.4 per 
cent credibility intervals for the (marginalized) joint probability distribution of two bias parameters. On the diagonal, we indicate the 
corresponding confidence levels for the one-parameter fits. The results for the low, medium and high redshiflt samples are plotted from 
left to right: bag clearly increases with redshift. Right panel: Same as left panel but using more principal components of the errors so as 
to account for 99 per cent of the variance (as opposed to 95 per cent in the left panel). Note that errorbars shrink and bp is shifted by 
almost 1(T with respect to the left panel. 



is such that at a given redshift we can only probe a fac- 
tor of ten in luminosity. This implies that any luminosity 
segregation has to be observed over this relatively small lu- 
minosity range. By slicing the different redshift bins into 
six sub-samples, our analysis tries to extract the maximum 
possible information from these limited data. However, to 
avoid to be dominated by measurement noise, we compare 
narrow luminosity intervals (at both the bright and faint 
ends) with larger control samples. This naturally leads to 
the selection of partially overlapping luminosity ranges. For 
this reason, our samples, despite probing different absolute 
magnitude intervals, are not independent from each other. 
Note the striking difference from previous studies of galaxy 
clustering (see e.g. Norberg et al. 2001, 2002) where one can 
perform several independent measurements of the cluster- 
ing amplitude as a function of luminosity and span a wide 
magnitude range. 

Given these limitations and based on theoretical mod- 
els of quasar clustering, should have we expected to detect 
some luminosity dependence in our dataset? As discussed 
in the introduction, models make very different predictions. 
From Fig. 3 in Lidz et al. (2006) we infer that the mod- 
els which assume a tight correlation between quasar lumi- 
nosity and the mass of the host-halos predict a difference 
A& = (tiB25 — &F25)o"8 ~ 1 for our high-redshift bin. On 
the other hand, models based on merger simulations predict 
A6 ~ 0.2. Similarly, the semi-analytic models by Kauffman 
& Haehnelt (2002) predict 0.3 < A6 < 0.7 depending on the 
characteristic quasar lifetime (this increases up to Ab ~ 1 
at 2: ~ 1.5). In all cases, the errorbars associated with our 
measurements are too large to lead to a statistically signif- 



icant detection of luminosity dependent clustering (in fact, 
our error on A6 is aAt — 2.3 and we measure Ab ~ 4.3 
which is only 1.4 ct and 1.8 a away from the two reference 
models). Assuming that Ab ~ 1, a sample which is nearly 
fifty times larger than ours (i.e. nearly all-sky) is needed 
to detect the luminosity dependence at the 3 a confidence 
level (assuming that the error scales with the sample size, 
n, as n~^^'^). Very likely a sample of this size will not be 
available in the next few years. Alternatively, one could try 
to beat the noise by reaching a fainter apparent magnitude 
limit. Fig. 3 in Lidz et al. (2006) shows that current models 
predict that, by reducing the median luminosity of our F25 
sample by a factor of 10, one would get 0.6 < A6 < 2. Note 
that a value of A6 ~ 2 could be ruled out at the 3 a confi- 
dence level by using a sample twelve times as large as ours. 
In conclusion, our analysis can only rule out (or detect) ex- 
treme luminosity dependent biasing and future surveys are 
required to discriminate among current models. 



5 REDSHIFT-DEPENDENT CLUSTERING 

In this section we focus on the redshift dependence of the 
quasar clustering amplitude. This is done in real space, using 
the projected correlation function. Our study is thus com- 
plementary to the redshift-space analysis by Groom et al. 
(2005) which might be affected by large-scale infall. In this 
section we also discuss the effect of using different methods 
for determining the quasar correlation length. 

To improve the analysis performed by PMN04 and bet- 
ter study the evolution with cosmic time, we split our quasar 
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Table 3. Best-fitting constant-bias (columns 4,5 and 6) and power-law models for six complementary redshift bins. The goodness of 

each fit is measured by the quantity x^jn/dof which gives the minimum value assumed by the statistic divided by the number of 
degrees of freedom. 

Zl 2h Zn, b n, xLiJdol ro 7 Xmin/dof f'^^^'**' Xmin/dof 

(/i-^Mpc) (h-iMpc) (h-^Mpc) 

1-0.37 3.47/3 4.3^42 l-TOlg:!^ 2.02/2 4.9l°;f 2.05/3 

toll S-etl'g 1-07/3 ^.itii l-84lg'e| 0.87/2 4:.9ti° 0.89/3 

tS'i 6.511J 1.05/3 6.lt^:? 1.93lg;i^ 1-44/2 5.61°;^ 1.54/3 

-0.39 6.711J 5.33/3 4.1^3.1 l-SgloJs 4.92/2 5.4±°:| 5.17/3 

iass 0.53/3 7.2^2.? 2.03lo:47 1.19/2 ^.3^1° 1.48/3 

-0.55 ll-Sii:? 2.76/3 8.8^49 l-82lo:4g 2.86/2 S.6t\i 2.87/3 



0.80 


1.06 


0.933 


1.57: 


1.06 


1.30 


1.185 


1.76: 


1.30 


1.51 


1.410 


2.13: 


1.51 


1.70 


1.602 


2.33: 


1.70 


1.89 


1.796 


3.02: 


1.89 


2.10 


1.987 


4.13" 



sample into six redshift bins. These are obtained by dividing 
into two equal sized parts (in terms of quasar number) the 
redshift bins used in the previous sections and in PNM04. 
We first compute the projected auto-correlation functions 
of all the quasars which lie in a given redshift bin. We then 
use the method described in tl3.2l to fit these data with the 
model given in equation The best-fitting models cor- 
respond to very good values of the statistic, indicating 
that our models accurately describe the spatial dependence 
of the measured correlations. Finally, we estimate the cor- 
relation length, ro, by determining the scale at which the 
best-fitting model has £,'^{r) — 1. For comparison, we also 
fit a power-law model ^'^(r) — (r/roY'- We distinguish two 
cases where we either allow both ro and 7 to vary or we fix 
7 = 1.8. Our results are reported in Table (together with 
the corresponding 68.3 per cent confidence levels) and plot- 
ted in Fig. 13 Note that the best-fitting values for the quasar 
correlation length derived from the power-law model are sys- 
tematically lower than the CDM ones (the two estimates are 
anyway consistent within Icr uncertainties). This is because, 
in a CDM model, the power-law index of the mass auto- 
correlation function is a function of the spatial separation. 
For instance, 7cff — dln^/dlnr ~ 1.3 at r = 5 Mpc, 
while 7off — 1.9 at r = 20 h"^ Mpc. Similarly, the projected 
correlation function, H(rx)/rx, has a slope varying between 
1.6 and 2.3 when 5 < r < 20 /i"^ Mpc. Therefore, a CDM 
model has a much lower correlation on large-scales with re- 
spect to a power-law model with the same ro (and 7 < 2). 
It is then clear from equation Q that a CDM model needs 
a higher overall normalization than a power-law model to fit 
the same projected correlation. 

The evolution of the quasar correlation length between 
0.8 < z < 2.1 can be approximately described by the linear 
relation ® 

ro = [7.3 + 5.2 (z - 1.5)] /i"^ Mpc (6) 

(this becomes ro — [6.1+3.6 (z — 1.5)] h^^ Mpc for the values 
of ro determined with the power-law fit to H(rx)/rx). For 
0.8 < z < 1.6, our results are in good agreement with Croom 
et al. (2005), but we find evidence for a stronger variation 

^ A nearly perfect interpolant between our measured points is 
ro = + ("yr)*] Mpc. This, however, ignores the statisti- 
cal uncertainties in the measure of ro (i.e. has «C 1)- 



at high-z (see Fig. I^J. ^ This is more evident when we use 
the CDM model to fit the data. The models by Kauffman 
& Haehnelt (2002) with a quasar lifetime of ~ 10^ yr match 
our data at z < 1.7. However, as already pointed out in 
PMN04, an increase in the lifetime by nearly an order of 
magnitude is required to reproduce the biased-CDM fits at 
2 ~ 2. 

Given that the mass autocorrelation function rapidly in- 
creases with cosmic time, the bias parameter of our quasars 
has to increase with z. The bottom panel of Fig.|5]shows the 
evolution of b. Note that, within the range of cosmological 
models allowed by observations, the quasar bias parameter 
scales linearly with erg. For this reason we decided to plot 
the product Ug = 6 cts which is independent on the assumed 
value for the amplitude of the linear power spectrum. Simi- 
larly to the correlation length, cTg seems to increase rapidly 
for z > 1.6 and the bias evolution is well approximated by 
the relation 

^<^s^l+(i^)\ (7) 

In CDM models, the clustering amplitude of massive dark- 
matter halos mainly depends on their mass. It is then inter- 
esting to see what halo masses correspond to the bias param- 
eters of the 2QZ quasars. In the bottom panel of Fig.|3 we 
plot the evolution of (7| for a number halo masses which is 
obtained using the model by Sheth & Tormen (1999). * The 
observed bias evolution is well reproduced by assuming that 
quasars are associated with halos of mass ~ 10" Mo- This 
is fully consistent with the more detailed analysis performed 
in PMN04. 



Our estimates of the uncertainties for ro and 7 are nearly a fac- 
tor of 2 larger than in Croom et al. (2005). This is due to a number 
of facts. On one hand, we find that blockwise bootstrap errors on 
the redshift-space correlation function are typically 30% smaller 
than on the projected correlation function. On the other hand, 
Croom et al. (2005) assume statistically independent Poisson er- 
ror bars for the redshift-space correlation function at separations 
smaller than 50 Mpc. In particular, neglecting correlations 
between points at different spatial separations results into smaller 
errors for the fitted parameters. 

* A concordance cosmology with eg = 0.9 is assumed here but 
results depend only slightly on erg. 
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Figure 5. Top: Redshift evolution of tlie quasar correlation 
length. Filled circles with errorbars indicate the best-fitting val- 
ues obtained using a biased CDM model. The short-dashed line 
is the analytical fit to the data given in equation JSJ. Filled trian- 
gles with errorbars refer to the best-fitting values obtained with 
a power-law model (where both 7 and rg are free to vary). These 
points have been slightly displaced to the right to improve read- 
ibility. The long-dashed line is the linear fit to the data given 
in the main text. Empty symbols are analogous to their filled 
counterparts but are obtained using larger redshift intervals as 
in PMN04. The shaded area marks the 68.3 confidence levels 
obtained by Croom et al. (2005). Bottom: As in the top panel 
but for the quasar bias parameter. Horizontal errorbars mark the 
16^^ and the 68*'' percentiles of the redshift distribution in a 
bin. The dashed line is the fitting function given in equation 171 . 
The dotted lines indicate the bias evolution for dark matter halos 
with masses lO^^ , lO^^'^^ , lO" , lO^^-^ and IO^Mq (from bottom 
to top). 

6 SUMMARY 

We have used a flux limited sample of ~ 14, 000 2QZ quasars 
brighter than Mbj — Slogj^Q h = —21.7 to study the quasar 
clustering properties in the redshift range 0.8 < z < 2.1. 
Our main results are summarized as follows. 

(i) Splitting the sample in three redshift intervals each 
divided into six luminosity ranges and combining infor- 
mation from the corresponding auto- and cross-correlation 
functions, we find some evidence for clustering segregation 
with luminosity. For redshifts z > 1.3, a frequentist model 
selection technique (the F-test) prefers a multi-parameter fit 
to the data at the 95 per cent confidence level 

(ii) A number of statistical tests based on information 
theory and Bayesian techniques show weak evidence for lu- 
minosity dependent clustering at high redshift {z > 1.3) and 
no evidence at low redshift {z < 1.3). These results some- 
what depend on the number of principal components used in 
the fitting procedure. Accounting for a larger fraction of the 
bootstrap variance increases the significance of the detection 
of clustering segregation with luminosity. 

(ill) Larger datasets, possibly with a deeper coverage, 
are needed to discriminate among current models of quasar 



formation and to pin-point the detailed quasar clustering 
trends as a function of luminosity at a given redshift. 

(iv) Splitting the sample into six complementary red- 
shift bins, we find strong evidence for an increase of the clus- 
tering amplitude with lookback time. We detect pure quasar- 
clustering evolution between Zoff = 0.93 and Zas = 1.99 at 
the 2.7a confidence level. A linear fit for the evolution of ro 
with redshift is given in equation (|SJ. 

(v) Accounting for the evolution of the mass density 
in a ACDM model, we find that the high-redshift quasars 
(zoff = 1.99) are ~ 2.6 times more biased than their low- 
redshift counterparts (zcff = 0.93). Evolution in b is detected 
at the 4.3(T confidence level. 

(vi) The clustering amplitude of optically selected 
quasars suggests that they are hosted by halos with mass 
M 10" M0 (see also PMN04). 
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APPENDIX A: MODEL SELECTION CRITERIA 

The F-test. If the uncertainties of the measurements are 
known to a good precision, model selection can be performed 
using a simple test. However, we cannot prove that our 
bootstrap errors accurately reproduce the true scatter in the 
data. An alternate test that does not require knowledge of 
the true standard deviation (up to a constant scaling factor) 
is to form the ratio 

^ Xn — fc Xn ^ ^ f A -l\ 

^ ~k ^ 

where A'' is the number of independent data values used 
in the fit (in our case A*' = A^pca)- The numerator in this 
equation is the difference between the calculated for n 
parameters and that for n — k parameters. Being a ratio 
of two x^-distributed variables, T follows an f -distribution 
with (fc, N — n) degrees of freedom. We can then estimate 
how significant is the addition of k parameters by integrating 
this distribution from to T. 

The Akaike information criterion. The question to find 
which model best approximate a given dataset can be ad- 
dressed in terms of information theory: the best model min- 
imizes the loss of information. The Akaike information cri- 
terion (AIC) evaluates models using the KuUback-Leibler 
information (Akaike 1973). In terms of the number of inde- 
pendent datapoints, n, and of model parameters, p, ^ 

AIC = xLn + 2p+^^^^^ (A2) 

n — p — i 

^ In our case p is the number of model parameters plus one to 
account for the estimation of the function. 
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where the second order term on the right is needed when 
the size of the dataset does not exceed the number of model 
parameters by a large factor (> 40). The model with the 
minimum AIC has to be considered the best. Note that the 
AIC penalizes for the addition of parameters and thus selects 
a model that fits well but has a minimum number of param- 
eters. The relative probability that a model is the correct 
solution is given by the Akaike weights 

-AAICi/2 

= —§ (A3) 

where the indices i and k run over the different models and 
A AIC = AICi — minfe(AICfc). In general, models receiv- 
ing AIC within 2 of the "best" deserve consideration, those 
within 3-7 have considerably less support and those above 
7 are basically rejected by the test (Burnham & Anderson 
1998). 

The Bayes factor. Bayes factors are the dominant method 
for Bayesian model testing (see Kass & Raftery 1994 for an 
extensive review). The Bayes factor, -B12, is the ratio be- 
tween the marginalized likelihoods (i.e. between the proba- 
bilities of the data given the models) for two different mod- 
els (1 and 2) and provides a scale of evidence in favour of a 
model versus another. The model with the highest marginal- 
ized likelihood is the best. Following an early suggestion by 
Jeffreys (1961), Kass & Raftery (1994) proposed the fol- 
lowing "rule of thumb" for interpreting the Bayes factors: 
1 < B12 < 3 provides weak evidence (barely worth men- 
tioning) for model 1, 3 < B12 < 12 provides substantial 
evidence for model 1, 12 < B12 < 150 provides strong evi- 
dence for model 1, B12 > 150 provides decisive evidence for 
model 1. Several numerical approaches have been proposed 
to compute Bayes factors based on MCMC sampling but 
most of these methods are subject to numerical or stability 
problems (see e.g. Han & Carling 2000 for a review). As a 
compromise between accuracy and coding complexity we es- 
timate the Bayes factors using the harmonic mean estimator 
(see e.g. Kass & Raftery 1994). We repeat the calculation 
using three different chains and we use the standard devi- 
ation among them as an estimate of the uncertainty of our 
marginalized likelihoods. 

The Bayesian information criterion. Unfortunately, while 
Bayes Factors are rather intuitive, as a matter of fact they 
are often quite difficult to calculate. This makes simpler (but 
approximate) estimates of the Bayes factors of great interest. 
Schwarz (1978) derived the Bayesian information criterion 
(BIC) as a large sample approximation to twice the loga- 
rithm of the Bayes factor. 

BICi - BIC2 ~ -21n(Bi2) , (A4) 

with 

BIC = xLn+plnn. (A5) 

Similarly to the AIC, the preferred model is the one with 
the lowest value of the criterion. The Kass-Raftery criterion 
for model selection is also applied to the BIC. 
The Deviance Information Criterion. Spiegelhalter et al. 
(2002) have recently proposed a Bayesian generalization of 
the AIC: the deviance information criterion (DIC). It is 
based on the posterior distibution of the log-likelihood or 
deviance {D = — 21og£ which, in our case, coincides with 



function): 

mC={D)+pD (A6) 

where (D) is the posterior expectation of the deviance while 
the effective number of parameters, po, is defined as the dif- 
ference between the posterior mean of the deviance and the 
deviance evaluated at the posterior mean of the parameters. 
Pd = (D) - D{{e)). The DIC can be re-written as 

mC = D{{e})+2pD (A7) 

which makes the analogy with the first-order AIC explicit. 
A good model corresponds to a low DIC and model-selection 
criteria developed for the AIC appear to work well also for 
the DIC (Spiegelhalter et al. 2002). Note, however, that, 
since pn < p, the DIC tends to be less conservative than 
the AIC in terms of model complexity. The main attraction 
of using this measure is that it is trivial to compute when 
performing MCMC on the models. 
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